Slow dynamics of water under pressure 
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We perform lengthy molecular dynamics simulations of the SPC/E model of water to investigate 
the dynamics under pressure at many temperatures and compare with experimental measurements. 
We calculate the isochrones of the diffusion constant D and observe power-law behavior of D on 
lowering temperature with an apparent singularity at a temperature Tc{P), as observed for water. 
Additional calculations show that the dynamics of the SPC/E model are consistent with slowing 
down due to the transient caging of molecules, as described by the mode-coupling theory (MCT). 
This supports the hypothesis that the apparent divergences of dynamic quantities along Tc{P) in 
water may be associated with "slowing down" as described by MCT. 



On supercooling water at atmospheric pressure, many 
thermodynamic and dynamic quantities show power-law 
growth (l). This power law behavior also appears under 
pressure, which allows measurement of the locus of ap- 
parent power-law singularities in water [Fig. |l|(a)]. The 
possible explanations of this behavior have generated a 
great deal of interest. In particular, three scenarios have 
been considered: (i) the existence of a spinodal bounding 
the stability of the liquid in the superheated, stretched, 
and supercooled states [Q; (ii) the existence of a hquid- 
liquid transition line between two liquid phases differing 
in density 1^-^; (iii) a singularity- free scenario in which 
the thermodynamic anomalies are related to the pres- 
ence of low-density and low-entropy structural hetero- 
geneities 1^ . Based on both experiments |p|,p|JlO|] and re- 
cent simulations ||ll[ , several authors have suggested that 
the power-law behavior of dynamic quantities might be 
explained by the transient caging of molecules by neigh- 
boring molecules, as described by the mode-coupling the- 
ory (MCT) ||l^, which we address here. This explanation 
would indicate that the dynamics of water are explain- 
able in the same framework developed for other fragile 
liquids |l3) , at least for temperatures above the homoge- 
neous nucleation temperature Tjj. Moreover, this expla- 
nation of the dynamic behavior on supercooling may be 
independent of the above scenarios suggested for thermo- 
dynamic behavior [Fig. |l|(a)]. 

Here we focus on the behavior of the diffusion constant 
D under pressure, which has been studied experimen- 
tally [^. We perform molecular dynamics simulations in 
the temperature range 210 K - 350 K for densities rang- 
ing from 0.95 g/cm"^ - 1.40 g/cm^ [Table |] using the ex- 
tended simple point charge potential (SPC/E) |l^. We 
select the SPC/E potential because it has been previ- 
ously shown to display power-law behavior of dynamic 
quantities, as observed in supercooled water at ambient 
pressure 15|. 

In Fig. ^, we compare the behavior of D under pressure 
at several temperatures for our simulations and the 
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FIG. 1. (a) Phase diagram of water. The extrapolated 
divergence of the isothermal compressibility (o) and the 
extrapolated divergence of D (filled o) M. The different 
locations of these divergences suggest that the phenomena 
may arise from different explanations, (b) Isochrones of D 
from simulation. The lines may be identified at follows: 
D = 10"^cmVs (o); D = lO"'^ '' cm^s (□); D = IQ-'^cmVs 
(o); D = 10"^cm^/s (A). The diffusion is also fit to 
D ~ {T/Tc - ly. The locus of Tc is indicated by (x). For 
reference, the (-|-) symbols indicate the locus of Tmd found 
in ref. ITgl. 
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FIG. 2. Diffusion constant D as a function of pressure for 
various temperatures from (a) our simulations and (b) NMR 
studies of water 

experiments of ref. The anomalous increase in D is 
qualitatively reproduced by SPC/E, but the quantitative 
increase of D is significantly larger than that observed ex- 
perimentally. This discrepancy may arise form the fact 
that the SPC/E potential is under- structured relative to 
water Jl9[ , so applying pressure allows for more bond 
breaking and thus greater diffusivity than observed ex- 
perimentally. We also find that the pressure where D 
begins to decrease with pressure - normal behavior for 
a liquid ~ is larger than that observed experimentally. 
This simple comparison of D leads us to expect that the 
qualitative dynamic features we observe in the SPC/E 
potential will aid in the understanding of the dynamics 
of water under pressure, but will likely not be quantita- 
tively accurate. 

We next determine the approximate form of the lines 
of constant D (isochrones) by interpolating our data over 
the region of the phase diagram studied [Fig. |(b)]. We 
note that the locus of points where the slope of the 
isochrones changes sign (i.e. the locus of points where 
D obtains a maximum value) is close to the Tmd lo- 
cus At each density studied, we fit Z? to a power 
law D ~ {T/Tc - 1)^. The shape of the locus of val- 
ues compares well with that observed experimentally , 
and changes slope at the same pressure [Figs. ||(a) and 
(b)]. We find the striking feature that 7 decreases under 
pressure for the SPC/E model, while 7 increases exper- 
imentally [Fig. ^. This disagreement underscores the 
need to improve the dynamic properties of water models, 
most of which already provide an adequate account of 
static properties [ pl[ . 

We next consider interpretation of our results using 
MCT, which has been used to quantitatively describe the 
weak supercooling regime - i.e., the temperature range 
where the characteristic times become three or 
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FIG. 3. Pressure dependence of the diffusivity exponent 7 
defined by 7? ~ {T/Tc — 1)'' • The symbols may be identified as 
follows: (o) 7 measured from simulation along isochores; (□) 
7 measured from simulation along isobaric paths, which are 
estimated from the isochoric data; (O) 7 measured along the 



; (A) experimental measurements 
It is clear from the available data 
that the SPC/E potential fails to reproduce the qualitative 
behavior of 7 under pressure in liquid water. 



-80 MPa isobar in ref. 
of 7 in water from ref. 



four orders of magnitude larger than those of the nor- 
mal liquid [ p^ . The region where experimental data 
are available in supercooled water is exactly the region 
where MCT holds. MCT provides a theoretical frame- 
work in which the slowing down of the dynamics arises 
from caging effects, related to the coupling between den- 
sity modes, mainly over length scales on the order of the 
nearest neighbors. In this respect, MCT does not require 
the presence of a thermodynamic instability to explain 
the power-law behavior of the characteristic times. 

MCT predicts power-law behavior of D, and also that 
the Fourier transform of the density-density correlation 
function F{q, t), typically referred to as the intermediate 
scattering function, decays via a two-step process. F{q, t) 
can be measured by neutron scattering experiments and 
is calculated via 



Fiq,t) 




-iq-[rfc(t)-rj(0)] 



(1) 



where S{q) is the structure factor [|2^. In the first relax- 
ation step, F{q, t) approaches a plateau value -F'piatcau('?); 
the decay from the plateau has the form -Fpiatcaul?) ~ 
F{q,t) ~ , where b is known as the von Schweidler ex- 
ponent. According to MCT, the value b is completely 
determined by the value of 7 [Q , so calculation of these 
exponents for SPC/E determines if MCT is consistent 
with our results. 
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FIG. 4. Fit of the stretched exponential of Eq. (^) for 
t > 2 ps at T = 210 K to both Fscif{q,t) (o) and F{q,t) 
(□) to obtain The horizontal line indicates the value pre- 
dicted by MCT for b using 7 values extrapolated from Fig. pi 
For P > 80 MPa, the relaxation of F(g, t) for q > 60 nm^ 
comes almost entirely from the first decay region, so the /3 
values obtained are not reliable in this range. Longer simula- 
tions, currently underway, will produce more reliable results 
in this region. 



The range of validity of the power-law t'' is strongly 
g-dependent making unambiguous calculation of b 
difficult. Fortunately, the same exponent b controls the 
long-time behavior of F{q, t) at large q. Indeed, MCT 
predicts that at long time, F(q,t) decays according to a 
Kohlrausch- Williams- Watts stretched exponential 



F{q,t) = A{q)e^p 



(2) 
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FIG. 5. The line shows the predicted relationship between 
b and 7 from MCT. The symbols show the calculated values 
for the SPC/E model; (o) from this work, (filled O) from 
ref. [0. 



is destroyed under pressure and that the water dynam- 
ics become closer to that of normal liquids, where core 
repulsion is the dominant mechanism. 

A significant result of our analysis is the demonstration 
that MCT is able to rationalize the dynamic behavior of 
the SPC/E model of water at all pressures. In doing so, 
MCT encompasses both the behavior at low pressures, 
where the mobility is essentially controlled by the pres- 
ence of strong energetic cages of hydrogen bonds, and 
at high pressures, where the dynamics are dominated by 
excluded volume effects. 

We wish to thank A. Rinaldi, S. Sastry, and A. Scala 
for their assistance. FWS is supported by an NSF grad- 
uate fellowship. The Center for Polymer Studies is sup- 
ported by NSF grant CH9728854 and British Petroleum. 



with limg^oo /3(g) = b ||2^. We show the g-dependence 
of f3 for each density studied at T = 210 K [Fig. |. We 
also calculate (3 for the "self-part" of F{q,t), denoted 
Fscu{q,t) (27|. In addition, we show the expected value 
of b according to MCT, using the values of 7 extrapolated 
from Fig. ||. The large-g limit of /? appears to approach 
the value predicted by MCT . Hence we conclude that 
the dynamic behavior of the SPC/E potential in the pres- 
sure range we study is consistent with slowing down as 
described by MCT [Fig. |]. We also note that on increas- 
ing pressure, the values of the exponents become closer 
to those for hard-sphere (7 = 2.58 and b = 0.545) and 
Lennard-Jones (7 = 2.37 and b = 0.617) systems [ p9[ . 
This confirms that the hydrogen-bond network 
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TABLE I. Summary of the state points simulated. We simulate 216 water molecules interacting via the SPC/E pair 
potential jl^. We simulate two independent systems at all temperatures (except 350 K), as the large correlation time makes 
time averaging more difficult. We equilibrate all simulated state points to a constant temperature by monitoring the pressure 
and internal energy. We control the temperature using the Berendsen method of rescaling the velocities with a thermostat 
time of 200 ps. The reaction-field technique with a cutoff of 0.79 nm accounts for the long-range Coulomb interactions [ p^ . 
The equations of motion evolve using the SHAKE algorithm |l| with a time step of 1 fs. Additional details can be found in 
ref. Systems are equilibrated for a time icq, followed by data collection runs for a time tdata- For all state points, the 

uncertainty in the potential energy U is less the 0.05 kj/mol. 
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